$title  DATAPREPARATION.GMS

* ###################################

* Per capita income, consumption patterns, and $CO_2$ emissions

* Journal of the Association of Environmental and Resource Economists

* Justin Caron and Thibault Fally

* April 2021

*###################################

* PREPARES IO COEFFICIENTS, TRADE SHARES, ETC, TO BE PASSED ON TO THE SIMULATION / STATISTICS FILES..


* load gravity data from stata or gams ?
* choices : stata, gams
$ if not set gravitydata $set gravitydata gams


$set SLASH \
$if %system.filesys% == UNIX $set SLASH /


$if not set datadir $set datadir "data%SLASH%"
$setglobal datadir %datadir%


* --------------------------------------------
* DEFINE CALIBRATION TYPE

* -- set of sectors - options : all, fdonly
$if not set sectorset $set sectorset all

* -- type of factor intensity - options : avg, cs - if using average, should change definiton of ENDOW
$if not set beta $set beta avg

* -- rescale ? - options : yes, yes
$if not set rescale $set rescale no

* declare set of demand estimates to use : theta4, notc, tc
$if not set demandest $set demandest tc

* -- declare dataset
$if not set ds $set ds gtap8gas


* define set of industries and countries here (overrides what is loaded from loaddata file)
Sets
r country  /
AUS,NZL,CHN,HKG,JPN,KOR,MNG,TWN,KHM,IDN,LAO,MYS,PHL,SGP,THA,VNM,BGD,IND,NPL,PAK,LKA,CAN,
USA,MEX,ARG,BOL,BRA,CHL,COL,ECU,PRY,PER,URY,VEN,CRI,GTM,HND,NIC,PAN,SLV,AUT,BEL,CYP,CZE,
DNK,EST,FIN,FRA,DEU,GRC,HUN,IRL,ITA,LVA,LTU,LUX,MLT,NLD,POL,PRT,SVK,SVN,ESP,SWE,GBR,CHE,
NOR,ALB,BGR,BLR,HRV,ROU,RUS,UKR,KAZ,KGZ,ARM,AZE,GEO,BHR,IRN,KWT,QAT,SAU,TUR,ARE,EGY,MAR,
TUN,CMR,CIV,GHA,NGA,SEN,ETH,KEN,MDG,MWI,MUS,MOZ,TZA,UGA,ZMB,ZWE,BWA,NAM,ZAF,ISR,OMN
/;

Sets
g   /
isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,omn,cmt,omt,vol,mil,pcr,
sgr,ofd,b_t,tex,wap,lea,lum,ppp,p_c,crp,nmm,i_s,nfm,fmp,mvh,otn,ele,ome,omf,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi, c,g, i
/;


Sets
i(g) industry  /
isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,omn,cmt,omt,vol,mil,pcr,
sgr,ofd,b_t,tex,wap,lea,lum,ppp,p_c,crp,nmm,i_s,nfm,fmp,mvh,otn,ele,ome,omf,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi
/;

$include loaddata_gtap8.gms



* ------------------------------------
* IMPORT GRAVITY AND DEMAND ESTIMATES

* gravity
parameter	coeffs
		ex,im,
		cst
		phiest;

* demand
parameter  lambdaest_;

parameter
fe           previously estimated demand fixed effects
sigma(i)        previously estimated sigmas
scale           scaling parameter for sigma
technology(i,r) export fixed effects from GE
tcost(i,r,s)    fitted transport costs from GE
theta           calibrated using estimates from the literature?
fittedexp	fitted expenditure
;


$gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_rall.gdx

$load  coeffs, ex, im, cst, fe,  lambdaest_=lambda, scale= sigma_scale fittedexp



* theta assumed to be 4
theta = 4;
sigma(i) = scale("non-homoth") * coeffs("sigma","coeff",i);


* ---------------------------------------------------------------
* ---------------------------------------------------------------


* ---------------------------------------------------------------
* -- DEFINE SET OF SECTORS AND REGIONS TO INCLUDE IN MODEL

set i_(i), r_(r), g_(g); i_(i) = no; r_(r) = no;

$if "%sectorset%"=="all"  i_(i) = yes;

*i_("dwe") = no;

$if "%sectorset%"=="fdonly"  i_(i)$sigma(i) = yes;

g_(i_) = yes; g_("c") = yes; g_("g") = yes; g_("i") = yes;

set c(g);
c("c") = yes; c("g") = yes; c("i") = yes;


set fds(i) set of sectors with final demand calc;

fds(i)$sigma(i) = yes;
*fds("dwe") = no;

* dropping the FDS structure
*fds(i_) = yes;
*sigma(i_)$(not sigma(i_)) = scale;

set countries(r) /

AUS,NZL,CHN,HKG,JPN,KOR,MNG,TWN,KHM,IDN,LAO,MYS,PHL,SGP,THA,VNM,BGD,IND,NPL,PAK,LKA,CAN,
USA,MEX,ARG,BOL,BRA,CHL,COL,ECU,PRY,PER,URY,VEN,CRI,GTM,HND,NIC,PAN,SLV,AUT,BEL,CYP,CZE,
DNK,EST,FIN,FRA,DEU,GRC,HUN,IRL,ITA,LVA,LTU,LUX,MLT,NLD,POL,PRT,SVK,SVN,ESP,SWE,GBR,CHE,
NOR,ALB,BGR,BLR,HRV,ROU,RUS,UKR,KAZ,KGZ,ARM,AZE,GEO,BHR,IRN,KWT,QAT,SAU,TUR,ARE,EGY,MAR,
TUN,CMR,CIV,GHA,NGA,SEN,ETH,KEN,MDG,MWI,MUS,MOZ,TZA,UGA,ZMB,ZWE,BWA,NAM,ZAF,ISR,OMN

/;


set     serv(i) service sectors / CMN, ISR,OBS, OFI,OSG,ROS,WTR,TRD,CNS, OTP, ATP, WTP,  ELY/
        tradables(i) the tradable sectors;

tradables(i) = yes;
tradables(serv) = no;
display tradables;

coeffs("manufsector","coeff", i_)$(not tradables(i_)) = 1;


r_(r) = yes;

alias(r_,s_); alias(i_,j_); alias(i_,k_);


* AGGREGATE TO SMALLER SET OF SECTORS FOR REPORTING

SET  ii_	 Sectoral classifications /

MAN	Manufactured and Processed Goods
ENE
SRV
TRN
AGR
*EIS

/;


$setglobal source implan

$eolcom !


SET	mapi(i,*)	Mapping from GTAP sectors to EPPA-10 sectors /
ATP.TRN,
B_T.AGR,
C_B.AGR,
CMN.MAN,
CMT.AGR,
CNS.MAN,
COA.ENE,
CRP.MAN,
OIL.ENE,
CTL.AGR,
*DWE.SRV,
ELE.MAN,
ELY.ENE,
FMP.MAN,
FRS.AGR,
FSH.AGR,
GAS.ENE,
PDR.AGR,
WHT.AGR,
I_S.MAN,
ISR.SRV,
LEA.MAN,
LUM.MAN,
 MIL.AGR,
MVH.MAN,
NFM.MAN,
NMM.MAN,
OAP.AGR,
OBS.SRV,
OCR.AGR,
OFD.AGR,
OFI.SRV,
OME.MAN,
OMF.MAN,
OMN.MAN,
OMT.AGR,
OSD.AGR,
OSG.SRV,
OTN.MAN,
OTP.TRN,
P_C.ENE,
PCR.AGR,
PFB.AGR,
PPP.MAN,
ROS.SRV,
SGR.AGR,
TEX.MAN,
TRD.SRV,
V_F.AGR,
VOL.AGR,
WAP.MAN,
WTP.TRN,
WTR.SRV
/;


* -----------------------------------------------------------------------
* RECONSTRUCT ESTIMATES FROM GRAVITY - FROM GRAVITYINGAMS

* EXPORTER FIXED EFFECTS - at the power 1/theta to have technology as TFP
technology(i_,r_) = exp((ex(i_,r_) + cst(i_))/theta);
* 

* TRANSPORT COSTS - at  power -1/theta to have transport cost as a cost
* = "fitted" transport costs from the Gravity equations
tcost(i_,r_,s_) = exp( -(
                -importdist(r_,s_,"ldist")      * coeffs("ldist","coeff",i_) 
                +importdist(r_,s_,"contig")        * coeffs("contig","coeff",i_) 
                +importdist(r_,s_,"comlang_off")   * coeffs("comlang_off","coeff",i_) 
                +importdist(r_,s_,"colony")        * coeffs("colony","coeff",i_) 
                +importdist(r_,s_,"homebias")      * coeffs("homebias","coeff",i_) )
                 / theta);


* -------------------------------------------------------------------
* COMPUTE INITIAL VALUES

set scn scenarios / initial, bmk, data, cf, "%chg"/;

* define all variables of interest
parameter	FinalD final demand,
		pcincome per-capita income (actually expenditures for now),
		lambda,
		X absorbtion,
		Trade includes domestic absorbtion,
		price,
		PriceIndex,
		Phi,
		fprice,
		SkillPremium,
		OutputQty
		factord factor demand,
		interd intermediate demand,
		output pre-tax output net of excluded intermediate inputs
		macro macro stats
		ByReg(*,scn,*) 
		ByRegSect
		Bilat
		ByRegIO;

Endow(f,r) =  pcexp(r)  * pop(r) * sum(i_,vfm(f,i_,r)) / sum((f.local,i_),vfm(f,i_,r));
* Such that the sum of all income from factors equal GDP (normalizing fprice to one)

* this is production as defined by value of inputs
byregsect("outputVal","data",i,r) = sum(f,vfm(f,i,r)*(1+rtf(f,i,r))) + sum(j_, (vdfm(j_,i,r)*(1+rtfd0(j_,i,r)) + vifm(j_,i,r)*(1+rtfi0(j_,i,r))));



set nooutput(i,r);

nooutput(i_,r_) = yes;
nooutput(i_,r_)$byregsect("outputVal","data",i_,r_) = no;

display nooutput;

byregsect("fdVal","data",i_,r_) = fd(i_,r_);
trade("data",i_,r_,s_) = btrade(i_,r_,s_);
trade("data",i_,r_,r_) = sum(g_,vdfm(i_,g_,r_) *(1+rtfd0(i_,g_,r_)));

byregsect("X","data",i_,r_) = sum(g_, vdfm(i_,g_,r_)+ vifm(i_,g_,r_));

ByReg("PCIncomeNom","data",r_) = pcexp(r_);


factord("data",f,i_,r_) = vfm(f,i_,r_)*(1+rtf(f,i_,r_));
interd("value","data",j_,i_,r_) =  vdfm(j_,i_,r_)*(1+rtfd0(j_,i_,r_)) + vifm(j_,i_,r_)*(1+rtfi0(j_,i_,r_));


* --------------------------------------------------
* CALCULATE INPUT SHARES - BASED ON SET OF INCLUDED GOODS ONLY

parameter      intersh		intermediate input coefficient (from i to j)
	       beta		country specific factor input coefficient ;

* AVERAGES 
* NOTE : DOING ONLY ACROSS SELECTED SET OF REGIONS CHANGES RESULTS QUITE A BIT

* factor shares :  
$if "%beta%"=="avg"  beta(f,i_,r)$sum(r.local, byregsect("outputVal","data",i_,r)) = sum(r.local,vfm(f,i_,r)*(1+rtf(f,i_,r)))/ sum(r.local, byregsect("outputVal","data",i_,r));
$if "%beta%"=="cs"   beta(f,i_,r)$byregsect("outputVal","data",i_,r) = vfm(f,i_,r)*(1+rtf(f,i_,r))/ byregsect("outputVal","data",i_,r) ;

* intermediate input shares :
$if "%beta%"=="avg"	intersh(i_,j,r)   = sum(r.local, vdfm(i_,j,r)*(1+rtfd0(i_,j,r)) + vifm(i_,j,r)*(1+rtfi0(i_,j,r))) / sum(r.local,byregsect("outputVal","data",j,r));
$if "%beta%"=="cs"	intersh(i_,j,r)$byregsect("outputVal","data",j,r) = (vdfm(i_,j,r)*(1+rtfd0(i_,j,r)) + vifm(i_,j,r)*(1+rtfi0(i_,j,r))) / byregsect("outputVal","data",j,r);



* check zero profit condition
parameter chkzp	 zero-profit check ;

chkzp(i_,r_) = round((sum(f, beta(f,i_,r_)) + sum(j_, intersh(j_,i_,r_)) - 1), 10);

display chkzp;

X("data - with average coeff",i_,r_) = fd(i_,r_) + sum((j_,s_), intersh(i_,j_,r_) * trade("data",j_,r_,s_));



* this computes the full MRIO coefficients matrix. 

* idea: need to do it just once, and be able to use them for any other embodied carbon / energy computation later on

parameter   trade_sh(i,r,s)  share of country r in country s's absorption (Import shares)
		beta_bilat_dir direct bilateral IO coefficients
	    beta_bilat_total total bilateral IO coefficients
	    beta_bilat_total_new total bilateral IO coefficients
	    beta_bilat_total_abs
	    beta_dev;


set	iter /0*10/;

*set	iter /0*1/;



* trade shares defined as: from r to s:
trade_sh(i,r,s)$sum(r.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(r.local , trade("data",i,r,s));


* ----------------------------------------------------------------------------
* -- CALCULATE MRIO CARBON and ENERGY INTENSITIES  - loop method
* ----------------------------------------------------------------------------

* inconsistency between model and data:
* data actually tells us which percentage of absorption is imported (i.e we know the vdfm / vifm split)
* here need bilateral  ghosh

parameter production  ;

* ignoring vdfm/vifm split which GTAP describes
* best to keep it model based

production(i,r) = sum(s,    trade_sh(i,r,s)* (sum(j,  vdfm(i,j,s) + vifm(i,j,s) )   + fd(i,s)));

parameter ghosh, ghosh_own ; 

ghosh_own(i,r,s)$production(i,r) =  fd(i,s) * trade_sh(i,r,s) / production(i,r);


* bilateral ghosh
*share of output of i going to j in r for fd in s 
*ghosh(i,j, r, s)$(production(j,r) and production(i,r)) =   (trade_sh(j,r,s) * fd(j,s) / production(j,r)) * (vdfm(i,j,r)) / production(i,r);

ghosh(i,j, r, s)$(production(j,r) and production(i,r)) =   (trade_sh(i,r,s) *  (vifm(i,j,s) + vdfm(i,j,s))) / production(i,r);
ghosh(i,"sum",r, "sum") = sum((j,s), ghosh(i,j, r, s)) + sum(s, ghosh_own(i,r,s));


* ---------------------------------------------------------------------
* LOADING NON-CO2 DATA


$gdxin '%datadir%GTP_NCO2_MMTCEQ_V8.gdx'

* source: https://www.gtap.agecon.purdue.edu/resources/res_display.asp?RecordID=4343


parameter 

NC_QOCO2E Non-CO2 emissions assoc. with output by industries-M. m.ton C-eq. 
NC_ENDWCO2E Non-CO2 emissions assoc. with endowment by industries-M. m.ton C-eq.
NC_TRADCO2E Non-CO2 emissions assoc. with input use by industries-M. m.ton C-eq.
NC_HHCO2E Non-CO2 emissions assoc. with input use by households-M. m.ton C-eq. ;


$load NC_QOCO2E       NC_ENDWCO2E      NC_TRADCO2E      NC_HHCO2E 

set gas  all gases /co2, N2O, CH4, FGAS/;
set nco2(gas)  non-co2 gases /N2O, CH4, FGAS /;


* -----------------------------------------------------------------------
* ENERGY SET DEFINITION

Set	ene_all(i) all energy goods / ely, coa, oil, gas, p_c /,
	ff(i)   Fossil   fuels       /p_c, oil, gas, coa/     ,
	pe(i)	primary energy        /coa,oil,gas/   
	se(i)	secondary energy        /coa,ely,gas, p_c/  
        e(i)    Energy inputs        /coa,p_c,oil,gas,ely/ ;



* ---------------------------------------------------------------------
* CO2 AND ENERGY COEFFICIENTS


parameter	
		bmkco2	sectoral co2 emissions at benchmark
		energy  energy use,
		energyint  energy intensity,
		co2lim	exogenous co2 limits
		co2emit co2 emissions by fuel and sector
		co2int  co2 intensity by sector
		co2coeff co2 emission coefficient by fuel and sector
		regco2	regional co2 emissions;


bmkco2(ff,i,r) = eco2(ff,i,r) ;
bmkco2(ff,g,r) = eco2(ff,g,r) ;
bmkco2(ff,"fd",r) = (eco2(ff,"c",r) + eco2(ff,"g",r)) ;


* NEED TO AVOID DOUBLE COUNTING HERE

*two ways to track energy:
*- primary energy: track back to coal, gas, cru only (set ELE and p_c to 0 as this is only transformation)
*- secondary energy: track back to gas, ele, p_c (set fosils to ele = 0 and crude to oil to 0)
*difference: the first one tracks differences in the fuel mix of electricity and energy intensity of oil refining


* energy is in Mtoe
energy("primary",i,g,r)$pe(i) = evdd(i,g,r) +  evdi(i,g,r);
energy("primary",i,"fd",r)$pe(i) = energy("primary",i,"c",r) + energy("primary",i,"g",r);

energy("secondary",i,g,r)$(not sameas(g,"ely") and not sameas(g,"p_c")) = evdd(i,g,r) +  evdi(i,g,r);
energy("secondary",i,"fd",r) = energy("secondary",i,"c",r) + energy("secondary",i,"g",r);


* DIRECT CO2 COEFFICIENTS   
* NOTE: this is eta in the text
co2coeff(ff,g_,r_)$(vdfm(ff,g_,r_) + vifm(ff,g_,r_)) =  bmkco2(ff,g_,r_) / (vdfm(ff,g_,r_) + vifm(ff,g_,r_)) ;
co2coeff(ff,"fd",r_)$byregsect("fdval","data",ff,r_)  =  (bmkco2(ff,"fd",r_) ) / byregsect("fdval","data",ff,r_) ; 
co2coeff(ff,"fd","all")$sum(r_,byregsect("fdval","data",ff,r_))  =  sum(r_, bmkco2(ff,"fd",r_) ) / sum(r_,byregsect("fdval","data",ff,r_)) ; 

* average co2 coefficient but not average price
co2coeff(i,"fd avg tech","all")$(sum(r_,energy("secondary",i,"fd",r_))) = (sum(r_,bmkco2(i,"fd",r_) )) / (sum(r_,energy("secondary",i,"fd",r_))); 

* divide by price (=multiply by Q divide by D)
co2coeff(i,"fd avg tech",r_)$byregsect("fdval","data",i,r_)  = co2coeff(i,"fd avg tech","all") /  byregsect("fdval","data",i,r_) *  energy("primary",i,"fd",r_);

* how much do energy prices actually vary?
parameter energyprices;

energyprices(i,"fd",r)$energy("secondary",i,"fd",r) = byregsect("fdval","data",i,r) /  energy("secondary",i,"fd",r);

energyprices(i,"fd with tax",r)$energy("secondary",i,"fd",r) = 
(
(vdfm(i,"c",r)*(1+rtfd0(i,"c",r)) + vifm(i,"c",r)*(1+rtfi0(i,"c",r)))
+ (vdfm(i,"g",r)*(1+rtfd0(i,"g",r)) + vifm(i,"g",r)*(1+rtfi0(i,"g",r)))
)
/  energy("secondary",i,"fd",r);



*X = co2  * P / D  (avg) = co2 / Q

*needs to multiply quantiy

*co2 = X * Q = X * D/P

*so, need to actually export X / P


* compute averages across sectors
* ACTUALLY ARE ALL VERY CLOSE TO AVERAGE:
co2coeff(ff,"all",r_)$sum(g_, vdfm(ff,g_,r_) + vifm(ff,g_,r_)) =  sum(g_, bmkco2(ff,g_,r_)) / sum(g_, vdfm(ff,g_,r_) + vifm(ff,g_,r_)) ;

* and now average across regions :
co2coeff(ff,"all","all")$sum((r_,g_), vdfm(ff,g_,r_) + vifm(ff,g_,r_)) =  sum((r_,g_), bmkco2(ff,g_,r_)) / sum((r_,g_), vdfm(ff,g_,r_) + vifm(ff,g_,r_)) ;

* ---------------------
* CO2 INTENSITY METRICS

* compute co2 intensity of output
co2int("co2","output",i_,r_)$byregsect("outputVal","data",i_,r_) =  sum(ff, bmkco2(ff,i_,r_)) / byregsect("outputVal","data",i_,r_);

* average cO2 intensity of output
co2int("co2","output",i_,"avg") =  sum((r_,ff), bmkco2(ff,i_,r_)) / sum(r_,byregsect("outputVal","data",i_,r_));




* NON-CO2 INTENSITY METRICS


* to fix: if keeping GTAP 7 values, need to think about what to do with missing regions
* ALSO when computing average MRIO intensities way at the bottom (only ok if get the GTAP8 values)

parameter nco2output  total nco2 associated with output per sector;

* intermediate inputs are fossil fuels and crp (chemicals) 
* factors are land and capital

* summing: intermediate, factor and output

nco2output(nco2,i,r) =   sum(j,NC_TRADCO2E(nco2, j,i,r))+  sum(f, NC_ENDWCO2E(nco2,f,i,r))
   + NC_QOCO2E(nco2, i,r);

nco2output(nco2,"all",r) = sum(i, nco2output(nco2,i,r));
nco2output(nco2,i,"all") = sum(r, nco2output(nco2,i,r));
nco2output(nco2,"all","all") = sum((i,r), nco2output(nco2,i,r));

*display nco2output;

* household emissions are so so small that I ignore them

co2int(nco2, "output",i,r)$byregsect("outputVal","data",i,r) = nco2output(nco2,i,r) / byregsect("outputVal","data",i,r);

co2int(nco2, "output",i,"avg")$sum(r$nco2output(nco2,i,r), byregsect("outputVal","data",i,r)) = sum(r, nco2output(nco2,i,r))  / sum(r$nco2output(nco2,i,r), byregsect("outputVal","data",i,r));


* summing over all gases
co2int("GHG", "output",i,r) = sum(gas, co2int(gas, "output",i,r));
co2int("GHG", "output",i,"avg") = (sum((nco2,r), nco2output(nco2,i,r)) + sum((r,ff)$(sum(nco2,nco2output(nco2,i,r))), bmkco2(ff,i,r)))  / sum(r$(sum(nco2,nco2output(nco2,i,r))), byregsect("outputVal","data",i,r));
		   

co2emit("data",ff,i_,r_) = bmkco2(ff,i_,r_);
co2emit("data",ff,"fd",r_) = bmkco2(ff,"fd",r_);

regco2("output","data", r_) = sum((ff,i_),co2emit("data",ff, i_,r_)) + sum((ff),co2emit("data",ff, "fd",r_)); 


* energy needs to be treated differently than Co2:
* primary energy intensity

* direct energy: = energy per $ of coal, gas, oil
* output: 

* direct energy demand: = direct demand for ele, coal, gas, p_c
* indirect energy demand = demand for all other goods time embodied energy

* Fossil   fuels       /p_c, oil, gas, coa/ 
* pe = primary energy / coal, gas, oil/


* ENERGY INTENSITIES:


* direct: energy per $
energyint("primary","direct",pe,r_)$sum(g_, (vdfm(pe,g_,r_) + vifm(pe,g_,r_))) = sum(g_, energy("primary",pe,g_,r_))
   / sum(g_, (vdfm(pe,g_,r_) + vifm(pe,g_,r_))) ;


energyint("primary","direct",pe,"all")$sum((r_,g_), (vdfm(pe,g_,r_) + vifm(pe,g_,r_))) = sum((r_,g_), energy("primary",pe,g_,r_))
   / sum((r_,g_), (vdfm(pe,g_,r_) + vifm(pe,g_,r_))) ;


* output, primary = how much primary energy demand per $ of output
energyint("primary","output",i_,r_)$byregsect("outputVal","data",i_,r_) =  sum(pe, energy("primary",pe,i_,r_)) 
    / byregsect("outputVal","data",i_,r_);


* secondary, direct:

energyint("secondary","direct",se,r_)$sum(g_, (vdfm(se,g_,r_) + vifm(se,g_,r_))) = 
   sum(g_, energy("secondary",se,g_,r_))
   / sum(g_, (vdfm(se,g_,r_) + vifm(se,g_,r_))) ;

energyint("secondary","direct",se,"all")$sum((r_,g_), (vdfm(se,g_,r_) + vifm(se,g_,r_))) = sum((r_,g_), energy("secondary",se,g_,r_))
   / sum((r_,g_), (vdfm(se,g_,r_) + vifm(se,g_,r_))) ;

* gas and coal are to be considered as secondary here too, so need to adjust:

energyint("secondary","direct","gas",r_)$sum(g_$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("gas",g_,r_) + vifm("gas",g_,r_))) 
 = 
sum(g_, energy("secondary","gas",g_,r_))
   / sum(g_$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("gas",g_,r_) + vifm("gas",g_,r_))) ;

energyint("secondary","direct","gas","all")$sum((r_,g_)$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("gas",g_,r_) + vifm("gas",g_,r_)))
 = sum((r_,g_), energy("secondary","gas",g_,r_))
   / sum((r_,g_)$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("gas",g_,r_) + vifm("gas",g_,r_))) ;


energyint("secondary","direct","coa",r_)$sum(g_$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("coa",g_,r_) + vifm("coa",g_,r_))) 
= sum(g_, energy("secondary","coa",g_,r_))
   / sum(g_$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("coa",g_,r_) + vifm("coa",g_,r_))) ;

energyint("secondary","direct","coa","all")$sum((r_,g_)$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("coa",g_,r_) + vifm("coa",g_,r_))) = 
sum((r_,g_), energy("secondary","coa",g_,r_))
   / sum((r_,g_)$(not sameas(g_,"ely") and not sameas(g_,"p_c")), (vdfm("coa",g_,r_) + vifm("coa",g_,r_))) ;


* direct, FD only

energyint("secondary","direct fd",se,r_)$sum(c, (vdfm(se,c,r_) + vifm(se,c,r_))) = 
   sum(c, energy("secondary",se,c,r_))
   / sum(c, (vdfm(se,c,r_) + vifm(se,c,r_))) ;

energyint("secondary","direct fd",se,"all")$sum((r_,c), (vdfm(se,c,r_) + vifm(se,c,r_))) = sum((r_,c), energy("secondary",se,c,r_))
   / sum((r_,c), (vdfm(se,c,r_) + vifm(se,c,r_))) ;

* gas and coal are to be considered as secondary here too, so need to adjust:

energyint("secondary","direct fd","gas",r_)$sum(c$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("gas",c,r_) + vifm("gas",c,r_))) 
 = 
sum(c, energy("secondary","gas",c,r_))
   / sum(c$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("gas",c,r_) + vifm("gas",c,r_))) ;

energyint("secondary","direct fd","gas","all")$sum((r_,c)$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("gas",c,r_) + vifm("gas",c,r_)))
 = sum((r_,c), energy("secondary","gas",c,r_))
   / sum((r_,c)$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("gas",c,r_) + vifm("gas",c,r_))) ;


energyint("secondary","direct fd","coa",r_)$sum(c$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("coa",c,r_) + vifm("coa",c,r_))) 
= sum(c, energy("secondary","coa",c,r_))
   / sum(c$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("coa",c,r_) + vifm("coa",c,r_))) ;

energyint("secondary","direct fd","coa","all")$sum((r_,c)$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("coa",c,r_) + vifm("coa",c,r_))) = 
sum((r_,c), energy("secondary","coa",c,r_))
   / sum((r_,c)$(not sameas(c,"ely") and not sameas(c,"p_c")), (vdfm("coa",c,r_) + vifm("coa",c,r_))) ;



* output, primary = how much primary energy demand per $ of output
energyint("secondary","output",i_,r_)$byregsect("outputVal","data",i_,r_) =  sum(se, energy("secondary",se,i_,r_)) 
    / byregsect("outputVal","data",i_,r_);


energyint("secondary","output",i_,"avg")$sum(r_, byregsect("outputVal","data",i_,r_)) =  sum((r_,se), energy("secondary",se,i_,r_)) 
    / sum(r_, byregsect("outputVal","data",i_,r_));



* -------------------------------------------------------------------------------
* IO COEFFICIENT COMPUTATIONS


* NOT NEEDED ? COMPUTE WORLD AVERAGE GHOSH MATRIX:

$ONTEXT

parameter ownfd  output of i that goes to i;
parameter ghosh; 

set cons(g) /"c", "g", "i"/;

ownfd(i_) = sum((cons,r_),vdfm(i_,cons,r_)+ vifm(i_,cons,r_)) ;

ghosh(i_,"fd") = ownfd(i_) / sum((r_,g_),  vifm(i_,g_,r_) + vdfm(i_,g_,r_));

ghosh(i_,j_) = sum(r_,vifm(i_,j_,r_) + vdfm(i_,j_,r_))  / sum((r_,g_),  vifm(i_,g_,r_) + vdfm(i_,g_,r_));

ghosh(i_, "sum") = sum(j_, ghosh(i_,j_)) + ghosh(i_,"fd");

display ghosh;

alias(j_,k_);

set iterg/ 1*5/;

parameter	ghosh_tot, ghosh_tot_new;

ghosh_tot(i_,j_) = ghosh(i_,j_);

loop(iterg,
	ghosh_tot_new(i_,j_) = (
			sum(r_,vifm(i_,j_,r_) + vdfm(i_,j_,r_)) + 
		           sum(k_,  ghosh_tot(k_,j_) * sum(r_,vifm(i_,k_,r_) + vdfm(i_,k_,r_)))) /
* div output i
	sum((r_,g_),  vifm(i_,g_,r_) + vdfm(i_,g_,r_));

	ghosh_tot(i,j) = ghosh_tot_new(i,j);

*= ( co2e(g,r) + sum(i, am(i,r)*vifm(i,g,r))
*		            + sum(i, ay(i,r)*vdfm(i,g,r)) )/vom(g,r);
*	dev(iter,r) = smax(g, abs(ay(g,r)-aynew(g,r)));
*	ay(g,r) = aynew(g,r);


);


execute_unload 'ghosh.gdx' ;

$OFFTEXT



* ---- COMPUTE INDIRECT FACTOR INTENSITITES.

$ONTEXT

parameter	Amatrix	A matrix expressed as proportion of production,
		avgAmatrix
		IA I-A matrix,
		IAinverse(i,j) temporary inverted I-A matrix,
		IAinverted(i,j,*)  inverted matrix; 

amatrix(i,j,r)$vom(j,r) = (vdfm(i,j,r)+vifm(i,j,r))/vom(j,r);

* create (I-A)^(-1) matrix
loop(r,
* create (I-A) matrix which will be inverted
IA(i,i)$vom(i,r) = (1-Amatrix(i,i,r)); 
IA(i,j)$(vom(i,r) and not sameas(i,j))= -Amatrix(i,j,r);
* invert this matrix:
execute_unload 'gdxforinverse.gdx' i,IA;
execute 'invert gdxforinverse.gdx i IA gdxfrominverse.gdx IAinverse';
execute_load 'gdxfrominverse.gdx' , IAinverse;
IAinverted(i,j,r) = IAinverse(i,j);
);

* this should be quite close to output 1 : it is..
co2int("co2","output 2",i_,r_)$byregsect("outputVal","data",i_,r_) =  sum(ff, co2int(ff,i_,r_)* amatrix(ff,i_,r_)) ;

$OFFTEXT



* COMPUTE MRIO EQUIVALENTS: 


* ----------------------------------------------------------------------------
* -- CALCULATE MRIO CARBON INTENSITIES  - USING MODEL -- NOT USED
* ---------------------------------------------------

* NOTE: SEEMS TO OVERLOAD SOLVER. Only works for autarky case
* Let's keep the code anyways

$ontext

variable beta_co2_total(j,i,r), dummy;

equation eq_totalbeta, edummy ;

* CAN DROP THE PRODUCTION TERMS
*eq_totalbeta(ene_all,i,r)$production(i,r)..  production(i,r) * beta_co2_total(ene_all,i,r) =e= production(i,r) * beta_co2_output(ene_all,i,r) + sum(j, beta_co2_total(ene_all,j,r) * intersh(j,i,r) *  production(i,r));

* autarky case use this
*eq_totalbeta(ene_all,i,r)..   beta_co2_total(ene_all,i,r) =e=  beta_co2_output(ene_all,i,r) + sum(j, beta_co2_total(ene_all,j,r) * intersh(j,i,r) );

* multi-regional case use this
eq_totalbeta(ene_all,i,r)..   beta_co2_total(ene_all,i,r) =e=  beta_co2_output(ene_all,i,r) + sum((s,j)$(trade_sh(j,s,r)> 1e-4), trade_sh(j,s,r) * beta_co2_total(ene_all,j,s) * intersh(j,i,r));

* a bit slow.. use CPLEX?  -- PATH crashed.  .. or revert back to previous method.
* even CPLEX doesn't work
* problem too large

edummy.. dummy =e= 0; 

*model invert_beta / eq_totalbeta.beta_co2_total/;

model invert_beta / eq_totalbeta, edummy/;


*beta_co2_total.FX(j,i,r) =1;

dummy.l = 0;


beta_co2_total.UP(ene_all,i,r) =100;
beta_co2_total.LO(ene_all,i,r) =beta_co2_output(ene_all,i,r);

beta_co2_total.L(ene_all,i,r) = beta_co2_output(ene_all,i,r);

invert_beta.iterlim = 10000;

invert_beta.reslim = 1000000;


*option lp=cplex;

option lp=conopt;

solve invert_beta minimizing dummy using lp;

$offtext


* ----------------------------------------------------------------------------
* --  MRIO coefficients  - loop method  -- SLOW, NOT USED
* ----------------------------------------------------------------------------


* -----------
* computing the total IO coefficients directly

$ontext

* NOTE: THIS SEEMS TO BE SLOW NO MATTER HOW I DO IT.. 
* solution: compute desired metrics directly without explicitely computing the total IO coefficients



* !! doing some data trimming here.
* this cuts it down from 25 million non-zeros to 1.5 million

beta_bilat_dir(i,j,r,s)$((intersh(i,j,s) * intersh(i,j,s) * trade_sh(i,r,s) > 1e-6)) =   intersh(i,j,s) * trade_sh(i,r,s);

beta_bilat_total(i,j,r,s) = beta_bilat_dir(i,j,r,s);





alias(i,k);
alias(s,ss);

loop(iter,


*intermd

*beta_bilat_total_abs(i,k,j,r,s) =   (sum((ss),  beta_bilat_dir(k,j,ss,s) *  beta_bilat_total(i,k,r,ss) ));



beta_bilat_total_new(i,j,r,s) = beta_bilat_dir(i,j,r,s) + (sum((k,ss),  beta_bilat_dir(k,j,ss,s) *  beta_bilat_total(i,k,r,ss) ));


*beta_bilat_total_new(i,j,r,s) = beta_bilat_dir(i,j,r,s) + (sum((k),  beta_bilat_total_abs(i,k,j,r,s) ));


beta_dev(iter,j,s)$sum((i,r)  , beta_bilat_total(i,j,r,s)) = sum((i,r)  , beta_bilat_total_new(i,j,r,s)) / sum((i,r)  , beta_bilat_total(i,j,r,s));

*beta_bilat_total_iter(iter,i,j,r,s) = beta_bilat_total(i,j,r,s);

beta_bilat_total(i,j,r,s)$(beta_bilat_total_new(i,j,r,s) >  1e-6) = beta_bilat_total_new(i,j,r,s);

);

$offtext



* ----
* computing country-specific 

parameter IncElast, IncElast_new, inc_elast_iter ;

IncElast("direct",i,r) = sigma(i)*(sum(i.local,  fd(i,r))) /
                                               (sum(i.local,  sigma(i) * fd(i,r) ));


IncElast("total",i,r) = sum(s, ghosh_own(i,r,s)  * IncElast("direct",i,s)) + sum((j,s), ghosh(i,j, r, s) * IncElast("direct",j,s)); 




loop(iter,


IncElast_new("total",i,r) = sum(s, ghosh_own(i,r,s)  * IncElast("direct",i,s)) + sum((j,s), ghosh(i,j, r, s) * IncElast("total",j,s)); 

IncElast("total",i,r) = IncElast_new("total",i,r);

inc_elast_iter("total",i,r,iter) = IncElast("total",i,r);

);


* weigh by country's share in final demand 
coeffs("total incelast - mean shares - fdweights","coeff", i) = sum((r), IncElast("total",i,r)* fd(i,r)) /  (sum((r),  fd(i,r))) ;

* weigh by country's share of total absorbption (production for now)
coeffs("total incelast - mean shares","coeff", i) = sum((r), IncElast("total",i,r)* production(i,r)) /  (sum((r),  production(i,r))) ;

coeffs("total incelast - mean shares","coeff", ii_) = sum(i_$mapi(i_,ii_), coeffs("total incelast - mean shares","coeff",i_) * sum(r_, production(i_,r_))) /
     sum(i_$mapi(i_,ii_),  sum(r_, production(i_,r_)));



* ----------------------------------------------------------
* COMPUTING MRIO CO2/ENERGY COEFFICIENT  -- USE THIS

* note: if changes are made here, they should also be made in the co2estimations.gms file



set		metric /mrio, fuel, energy/;
parameters	ay(*,g,r)		MRIO carbon content of output
		ay2(g,r)	to compute SECONDARY energy content
		af(g,r)		Direct carbon content (fuel only)
		ad(g,r)		Direct carbon content (fuel + electricity)
		am(*,i,r)		Carbon content of imports
		am2
		at(j)		Carbon content of trade,
		at2
		atf(j)		Carbon content of trade (fuel only),
		atd(j)		Carbon content of trade (fuel + electricity)
		co2e(*,*,r)	Total carbon emissions by sector and region;



* adding non-co2 gases to this
co2e("co2",j,r) = sum(i, eco2(i,j,r));
*co2e("co2",g,r) = sum(i, eco2(i,g,r));

* non-co2 gases
co2e(nco2,j,r) = nco2output(nco2,j,r) ;

* all GHG
co2e("GHG",j,r) = co2e("co2",j,r) + sum(nco2, co2e(nco2,j,r));

* OVER WHICH GASES TO INVERT (takes a bit of time)
* do them all
*set gas_ /co2, nh4, n20, fgas, ghg/;

* only co2 and total GHG
set gas_ /co2, ghg/;



*co2e("co2","total",r) = sum(g,co2e("co2",g,r));

ay(gas,i,r)$production(i,r) = co2e(gas,i,r) / production(i,r);

ay2(j,r)$production(j,r) = sum(i,energy("secondary",i,j,r)) / production(j,r);

parameter	dev		Maximum deviation in embodied carbon
		aynew		New estimate of embodied carbon,
		ay2new
		ab(i,r,s,*,metric)	Carbon content of bilateral trade flow (gross and net);

*PARAMETER AY3, AY3NEW;

loop(iter,

* NOTE: ignoring the VIFM/VDFM split from GTAP and assuming identical trade shares everywhere
	aynew(gas_,i,r)$production(i,r) =  co2e(gas_,i,r) / production(i,r) + sum((j,s), trade_sh(j,s,r) *  intersh(j,i,r) * ay(gas_,j,s))	;

	dev(gas_,iter,r) = smax(g, abs(ay(gas_,g,r)-aynew(gas_,g,r)));
	ay(gas_,g,r) = aynew(gas_,g,r);

	ay2new(i,r)$production(i,r) =  sum(j,energy("secondary",j,i,r)) / production(i,r) + sum((j,s), trade_sh(j,s,r) *  intersh(j,i,r) * ay2(j,s))	;
	ay2(g,r) = ay2new(g,r);

);

* import intensity
am(gas_,i,r)$sum(s$(not sameas(r,s)),  trade_sh(i,s,r)) =  sum(s$(not sameas(r,s)),  trade_sh(i,s,r) *ay(gas_,i,s))/sum(s$(not sameas(r,s)),  trade_sh(i,s,r));


am2(i,r)$sum(s$(not sameas(r,s)),  trade_sh(i,s,r)) =  sum(s$(not sameas(r,s)),  trade_sh(i,s,r) *ay2(i,s))/sum(s$(not sameas(r,s)),  trade_sh(i,s,r));



* compute full MRIO co2 coeffs:
co2int(GAS_,"MRIO output",g,r_) = ay(gas_,g,r_);
co2int(gas_,"MRIO cons indirect",i,r_) = sum(s,  trade_sh(i,s,r_) *ay(gas_,i,s));
co2int("co2","MRIO cons total",i,r_) = co2int("co2","MRIO cons indirect",i,r_) + co2coeff(i,"fd",r_);
co2int(gas_,"MRIO imports",i,r_) = am(gas_,i,r_);


* output weighted by output
co2int(gas_,"MRIO output",i,"avg") = sum(r, production(i,r)* co2int(gas_,"MRIO output",i,r))/ sum(r,production(i,r));

* consumptin weighted by fd
co2int(gas_,"MRIO cons indirect",i,"avg") = sum(r, fd(i,r)* co2int(gas_,"MRIO cons indirect",i,r))/ sum(r,fd(i,r));
co2int(gas_,"MRIO cons total",i,"avg") = sum(r, fd(i,r)* co2int(gas_,"MRIO cons total",i,r))/ sum(r,fd(i,r));


* secondary energy

energyint("secondary","MRIO output",i,r) = ay2(i,r);
energyint("secondary","MRIO output",i,"avg")  = sum(r, production(i,r)* energyint("secondary","MRIO output",i,r))/ sum(r,production(i,r));

energyint("secondary","MRIO cons indirect",i,r) = sum(s,  trade_sh(i,s,r) *ay2(i,s));
energyint("secondary","MRIO cons indirect",i,"avg") = sum(r, fd(i,r)* energyint("secondary","MRIO cons indirect",i,r))/ sum(r,fd(i,r));

energyint("secondary","MRIO cons total",i,r) = sum(s,  trade_sh(i,s,r) *ay2(i,s))+ energyint("secondary","direct fd",i,r)  ;
energyint("secondary","MRIO cons total",i,"avg") = sum(r, fd(i,r)* energyint("secondary","MRIO cons total",i,r))/ sum(r,fd(i,r));


* ------------------
* FOR CORRELATION GRAPHS
    
* store in coeffs parameter to make the graphs
coeffs("CO2 int - fd direct","coeff", ff) = co2coeff(ff,"fd","all");

* GHG direct = co2 pretty much

coeffs("Co2 int - output","coeff", i) = co2int("co2","output",i,"avg") ;

coeffs("GHG int - output","coeff", i) = co2int("ghg","output",i,"avg") ;

coeffs("Sec en int - output","coeff", i) = energyint("secondary","output",i,"avg") ;


* ---------

* direct coeff  (correlate with TOTAL income elasticity)
coeffs("Co2 int - direct","coeff", i) = co2int("co2","output",i,"avg") + co2coeff(i,"fd","all");

* here, ignore non-co2 fd emissions (very very small)
coeffs("GHG int - direct","coeff", i) = co2int("ghg","output",i,"avg") + co2coeff(i,"fd","all");

* secondary direct  = output + fd direct 
* TO DO! USE ACTUAL FD.. NOT DIRECT (SHOULD BE ALMOST IDENTICAL THOUGH)
coeffs("Sec en int - direct","coeff", i) = energyint("secondary","output",i,"avg")  +energyint("secondary","direct",i,"all");


* weighted by sector's share of production emissions (incl. fd)

coeffs("sharedirectco2","coeff", i) =  (sum((r,ff),   bmkco2(ff,i,r)) + sum(r, bmkco2(i,"fd",r))) / sum(i.local, (sum((r,ff),   bmkco2(ff,i,r)) + sum(r, bmkco2(i,"fd",r))));
	    			    	  
coeffs("sharedirectGHG","coeff", i) =  ( sum((r,nco2),   nco2output(nco2,i,r)) + sum((r,ff), bmkco2(ff,i,r)) + sum(r, bmkco2(i,"fd",r))) /
		 sum(i.local,( sum((r,nco2),   nco2output(nco2,i,r)) + sum((r,ff), bmkco2(ff,i,r)) + sum(r, bmkco2(i,"fd",r))));




* ------------
* total coeff  (correlate with DIRECT income elasticity)
coeffs("Co2 int - total","coeff", i) = co2int("co2","MRIO output",i,"avg") + coeffs("CO2 int - fd direct","coeff", i);
					    
coeffs("GHG int - total","coeff", i) = co2int("ghg","MRIO output",i,"avg") + coeffs("CO2 int - fd direct","coeff", i);

coeffs("Sec en int - total","coeff", i) = energyint("secondary","MRIO output",i,"avg") + energyint("secondary","direct",i,"all");


* aggregate:
*coeffs("Co2 int - direct","coeff"	, ii_) = sum(i_$mapi(i_,ii_), coeffs("Co2 int - direct","coeff", i_) * sum(r_, fd(i_,r_))) /
*     sum(i_$mapi(i_,ii_),  sum(r_, fd(i_,r_)));

* better (note, does not include final co2 in cons)
coeffs("Co2 int - direct","coeff", ii_)  =  sum((r_,ff), sum(i_$mapi(i_,ii_),    bmkco2(ff,i_,r_))          ) / sum(r_,  sum(i_$mapi(i_,ii_),  byregsect("outputVal","data",i_,r_)));


coeffs("Co2 int - total","coeff", ii_) = sum(i_$mapi(i_,ii_), coeffs("Co2 int - total","coeff", i_) * sum(r_, fd(i_,r_))) /
     sum(i_$mapi(i_,ii_),  sum(r_, fd(i_,r_)));


*  shares in total consumption estimates for weights

* same as sh_co2_totalmrio("co2", i,r)  in next file
coeffs("sharetotalco2cons","coeff", i) =  sum(r, fd(i,r)* ( co2coeff(i,"fd",r) + co2int("co2","MRIO cons indirect",i,r))) / sum((r,i.local), fd(i,r)* ( co2coeff(i,"fd",r) + co2int("co2","MRIO cons indirect",i,r)) ) ;

* aggregate:
coeffs("sharetotalco2cons","coeff", ii_) = sum(i_$mapi(i_,ii_), sum(r, fd(i_,r)* ( co2coeff(i_,"fd",r) + co2int("co2","MRIO cons indirect",i_,r)))) / sum((r,i.local), fd(i,r)* ( co2coeff(i,"fd",r) + co2int("co2","MRIO cons indirect",i,r)) ) ;



	    			    
coeffs("sharetotalGHGcons","coeff", i) = sum(r, fd(i,r)* ( co2coeff(i,"fd",r) + co2int("ghg","MRIO cons indirect",i,r))) / sum((r,i.local), fd(i,r)* ( co2coeff(i,"fd",r) + co2int("ghg","MRIO cons indirect",i,r)) ) ;

coeffs("sharetotalsecencons","coeff", i) = sum(r, fd(i,r)* ( energyint("secondary","MRIO output",i,r) + energyint("secondary","direct",i,r))) 
		 / sum((r,i.local), fd(i,r)* ( energyint("secondary","MRIO output",i,r) + energyint("secondary","direct",i,r)) ) ;


* preparing other shares:

* share of total fd
coeffs("sharetotalfd","coeff", i) =  sum(r, fd(i,r)) / sum((r,i.local), fd(i,r) ) ;

* average share of fd
coeffs("sharefd avg","coeff", i) =   sum(r, fd(i,r) / sum((i.local), fd(i,r) ) ) / card(r);


execute_unload 'estimates%SLASH%co2_data_%demandest%_%beta%.gdx',  beta, intersh, co2int, energyint,  co2coeff, trade, fd, pcexp, incelast, bmkco2, energy,
 production, i, r, r_, i_, trade_sh, pop, nco2output, coeffs, energyprices, rtfd, rtfi, fd_disagg, endow, vom, vdfm, vifm, fd, vfm;





